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Abstract 

Traction forces exerted by adherent cells on their microenvironment can mediate many critical cellular functions. Accurate 
quantification of these forces is essential for mechanistic understanding of mechanotransduction. However, most existing 
methods of quantifying cellular forces are limited to single cells in isolation, whereas most physiological processes are 
inherently multi-cellular in nature where cell-cell and cell-microenvironment interactions determine the emergent 
properties of cell clusters. In the present study, a robust finite-element-method-based cell traction force microscopy 
technique is developed to estimate the traction forces produced by multiple isolated cells as well as cell clusters on soft 
substrates. The method accounts for the finite thickness of the substrate. Hence, cell cluster size can be larger than substrate 
thickness. The method allows computing the traction field from the substrate displacements within the cells' and clusters' 
boundaries. The displacement data outside these boundaries are not necessary. The utility of the method is demonstrated 
by computing the traction generated by multiple monkey kidney fibroblasts (MKF) and human colon cancerous {HCT-8) 
cells in close proximity, as well as by large clusters. It is found that cells act as individual contractile groups within clusters for 
generating traction. There may be multiple of such groups in the cluster, or the entire cluster may behave a single group. 
Individual cells do not form dipoles, but serve as a conduit of force (transmission lines) over long distances in the cluster. 
The cell-cell force can be either tensile or compressive depending on the cell-microenvironment interactions. 
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Introduction 

Recent research has demonstrated that cells communicate with 
each other as well as with their microenvironments through 
mechanical signaling [1,2,3,4,5,6], in addition to biochemical ones 
[7,8,9,10,11,12,13,14]. Many physiological processes, including 
cell adhesion [15,16,17], cytoskeleton polarity [13,18], cell 
proliferation [19,20], cell differentiation [12,21,22], embryogenesis 
[23,24], cancer metastasis [7,25], and wound-healing [26,27], can 
be significantly influenced by the transmission and sensation of 
physical forces between the cells and their microenvironments. For 
example, exposure of HCT-8 human colon cancer ceUs to soft 
substrates results in a profound stable cell state transition from an 
epithelial phenotype to a metastasis-like phenotype (MLP) 
[7,8,28,29,30,31]. Adherent ceUs actively sense the local anisotro- 
py of their microenvironment [2,18,32,33] as weU as the forces 
applied by neighboring cells [1,4,11,34,35], foUowed by polariza- 
tion of stress-fibers and synergetic cell functions. Hence, accurate 
estimation of the traction forces exerted by the cells on their 
substrates under various physiological conditions can provide 
important insight on many fundamental questions regarding the 
mechanical interactions between various cell types and their 



microenvironment [36,37,38]. Over the past few decades, several 
seminal techniques to assess the cellular traction forces have been 
developed (see reviews [14,39,40,41,42,43,44]). However, most of 
them are limited to computation of traction forces exerted by 
single, isolated cells. 

Efforts at visualizing cellular traction forces may be traced back 
to 1980s when Harris et al. used thin polymeric silicone substrates 
for cell culture, and observed the wrinkling phenomena caused by 
the traction of migrating cells [45]. However, quantitative 
estimation of the traction from the wrinkling of silicone substrates 
is challenging due to the inherent non-linearity of the problem. 
From 1995 on, Lee, Jacobsen and Dembo et al., as weU as other 
groups, developed several traction force microscopy techniques 
(TFM) to quantify the cellular traction produced by migrating or 
stationary cells on soft substrates [46,47,48,49,50,51,52,53,54]. 
TFM computes the cell traction forces from the deformation of a 
soft substrate with known elastic properties, such as polyacryl- 
amide (PA) gel, on which ceUs are cultured. The deformation is 
measured from the displacements of micro-fluorescent markers 
embedded in the substrate. The motion is measured from two 
images. First image is taken with the cells adhered to the substrate. 
Here, the cells have generated traction force on the substrate, and 
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Author Summary 

Adherent cells sense, transduce and respond to their 
microenvironment by generating traction forces on their 
surroundings. To accurately understand these mechano- 
transduction processes, it is critical to have a robust and 
reliable method for traction force visualization and 
quantification. However, most cell traction force micros- 
copy methods are limited to only single cell traction force 
analysis. Considering that most physiological processes are 
essentially collective multi-cellular events, there is a need 
for traction force microscopy methods capable of analyz- 
ing traction forces resulting from multiple cells. We have 
developed a novel and robust multi-cellular traction force 
microscopy method for computing cell traction on soft 
substrates, and applied it to compute traction field 
generated by both multiple cells and cell clusters. We 
verified the accuracy, robustness, and efficiency of the 
method by theoretical, numerical and experimental 
approaches. Our method provides a powerful toolset to 
pursue the mechanistic understanding of collective bio- 
logical activities, such as cancer metastasis and neuromus- 
cular interactions. 



the image gives the deformed configuration of the soft substrate. 
Then cells are removed from the substrate through trypsinzation, 
and a second image is taken. Subsequently, the substrate is 
relieved of cell traction, and the image shows the un-deformed 
configuration of the substrate. A comparison of the two images 
gives the displacement field of the substrate's top surface due to 
cell tractions. Digital image correlation method (DICM) is used to 
quantify the displacement field. The traction field is estimated 
from the displacement field. Several methods have been proposed 
for force estimation ranging from analytical methods, i.e. the 
Boussinesq formulation (either using Bayesian likelihood regular- 
ization method [51,55] or Fourier transformed approach [49]), to 
computational methods like finite element analysis (FEA) [56]. 
The Boussinesq formulation approach, which assumes the 
substrate as a semi- infinite elastic half space [57], was first adopted 
by Dembo and Wang, et al, to compute the traction forces from 
the displacement fields followed by regularization [51,55,58,59]. 
Since the Boussinesq formulation involves solving an inverse 
problem, the solution demands computational regularization 
schemes to predict the approximate traction solutions. Important- 
ly, Butler, Trepat and Fredberg, et al. [49,60,61,62,63] made 
significant progress in mitigating some pitfalls of the regularization 
scheme by solving the Boussinesq equation using Fourier 
transform. Later Schwarz et al. introduced a new method to 
compute traction forces only at the focal adhesion site of the cell by 
assuming that the cell force transfer occurs only through these 
sites [50]. Some novel platforms, such as the photobleaching- 
activated monolayer with adhesive micro-patterns developed by 
Scrimgeour et al. [64] and the elastic substrates with micro-contact 
printing demonstrated by Strieker et al. [65], were also used to 
characterize the cell traction force. Furthermore, a FEA-based 
technique was also developed by Yang et al. to greatly improve the 
accuracy of traction force calculations [56] . The FEA method no 
longer depends on the Boussinesq formulation and thus is not 
limited by the semi-infinite elastic half space assumption [66,67] . 
Recently additional contribution has been made in traction force 
computation in three dimensions [19,68,69,70,71,72,73]. 3D 
TFM techniques compute the 3D traction force fields from the 
cell induced 3D displacement and strain fields obtained using laser 
scanning coiifocal microscopy (LSCM) and digital volume 



correlation (DVC). However, it is challenging to obtain the Z- 
dimension displacement field and the technique can only be 
applied to single cell cases, rather than multiple cells or cell 
clusters. 

The above studies focused on traction force computation for 
single cells far from their neighbors, i.e. cells that do not interact 
mechanically with each other. However, live cells do interact with 
their neighbors chemo-mechanically and form cell clusters 
[7,29,37,74,75,76]. In this paper we present a novel finite- 
element-based TFM technique to compute the traction fields 
generated by multiple cells and clusters. We first present a 
theoretical proof showing that the 3D traction field computed 
from prescribed displacement field of the substrate is unique. We 
verify the uniqueness by considering a 2-ceU case. We test the 
accuracy of the computational technique by applying a known 
force on PA gel substrate using a micro-needle, and by comparing 
the experimental force with the computed one. Finally, we 
compute the traction fields generated by multiple cancerous and 
fibroblast cell clusters, and reveal that cells might be under 
compression in such 2D clusters. We believe that the present 
technique may enable better examination and understanding of a 
variety of biological phenomena involving homotypic and 
heterotypic cells and cell cluster interactions [77,78,79]. 

Results 

Uniqueness of traction field computed from 
displacement field in 3D linear elastic solids 

Consider a 3D linear elastic solid with volume V in static 
equilibrium. Its boundary, S, consists of >?„ and {S= S„+S^ 
where displacements uf and traction tf are prescribed respective- 

ly- 

Proposition: Given displacement field uf at S„ and traction tf at S^,, the 
corresponding traction f,- at is unique. (Note: indices i, j = 1,2,3 
correspond to x,y,z Cartesian coordinates respectively; all equa- 
tions follow standard tensor notation and summation convention). 
Supporting material Text S 1 presents the proof of the proposition. 

Simple ID examples of uniqueness 

Displacement boundary condition. To gain an intuitive 
insight on the uniqueness of traction solution, we present a simple 
ID model. The stiffness and compliance matrices of a uniform 
elastic bar have been derived (see Supplementary Materials Text 
S2). In Fig. la, a uniform bar is subjeted to three concentrated 
force Fi, F^, F^j at nodes 1, 2, 3 with corresponding displacements 
Uj, U2, and linear stiffness kj, k2, k^j respectively. For simphcity, 
let k] = 1, k2 = a, = h, then the displacement-force relationship is 
given by: 
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Thus, if displacements [uj, u^, M ;) are given, nodal forces [Fj, F2, 
F-i) can be obtained from Eqn (1). For a given displacement {ui, U2, 
U3j, there is only one possible value of (i^;, F2, F-^), i.e., solution for 
the nodal forces is unique, since Young modulus is positive, 
compliance and stiffness matrices are positive definite and hence 
are non-singular and invertible. In other words, there is always a 
unique relationship between displacement and forces on the 
nodes. 
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Figure 1. Modeling of cell contraction on 1D elastic substrate with mixed boundary conditions, (a) An elastic bar is discretized with 3 
nodes with concentrated forces applying on each node along with their respective displacements. Note that this general loading is used for deriving 
stiffness matrix which uniquely relates nodal forces to the nodal displacement subject to any boundary condition, (b) A cell applies contractile forces 
on nodes 1 and 2 (i.e with known, measured displacements) while node 3 is free. This set of inputs constitutes a Mixed Boundary Condition, in that a 
combination of nodal displacements and forces are given (i/i =0,»2 = — 1,^3 =0) and their respective unknowns {Fi,F2,u^) are computed by the 
model. 

doi:1 0.1 371/journal.pcbi.1 003631 .gOOl 



Mixed boundary value problem. In Eq (1), given all 
displacement data at the nodes, the force vectors can be directiy 
calculated, and vice versa. However, there are cases where a 
combination of displacement and forces on boundaries are given, 
called mixed boundary value problems (MBVP). For example, suppose a 
cell is adhered at nodes 1 and 2 and creates contractile forces Fj 
and F^. The corresponding displacements Uj and are measured 
but there is no anchorage at node 3 (i.e zero traction) and hence 
F-j-0 (Fig. lb). Given: uj = 0, U2 = —1, F:j = 0, let us find 
contractile forces Fj, F2 and displacement My. Our unknowns are 
a combination of forces and displacements. Applying the given 

boundary conditions into Eqn (1) and solving for the unknowns, 
a a 

we obtam: Fi = ,i^2 = — ,^3 = — 1 . 

I— a I— a 

Note that FjT^O although Uj = 0 and similarly u-,-¥= 0 while Fy = 0. 
Therefore, zero displacement does not necessarily result a zero 
force (traction) at a node and vice versa. The displacement at zero- 
traction node 3 is due to the displacements at other nodes. This 
example illustrates the counter intuitive possibility of non-zero 
displacements at points on the body where traction is zero. Also 
note that, Fi + = 0 which confirms the traction field under the 
cell is self-equilibrated. 

Finite element approach for solving cell traction on 2D 
substrates 

We illustrate our computational scheme as follows. Consider 
two separate cells on a soft elastic substrate. The substrate is 
adhered to a rigid surface (such as glass) at the bottom. The lateral 
boundary of the substrate is far from the cells. In the finite element 
scheme, the substrate is modeled as a rectangular pyramidal soUd 
body. It is discretized as a collection of small cubes with common 



nodes. We need to prescribe three boundary conditions, namely 
any combination of forces {F„ Fy, F^ and displacements {u„ Uy, mJ, 
at each of the surface nodes. For example, (i^, Uy, can be a 
boundary condition at a surface node. To ensure that the body is 
at rest (no rigid body translation or rotation), at least two of the 
nodes are prescribed with u„ = Uy, =u^ = 0. Given the boundary 
conditions, finite element scheme calculates the deformation of the 
solid body such that the total energy is minimized. Thus the 
displacements at each node within the body, and at the surface 
nodes where forces are prescribed are evaluated. This leads to the 
evaluation of strains and stresses using the elastic properties of the 
solid (Young's modulus and Poisson's ratio for the isotropic gel). 
Surface traction is calculated from the stress near the surface and 
normal vector to the surface (?,■ = aynj), as shown in Supplemen- 
tary Materials S 1 . Surface nodal forces are calculated from an area 
integral of traction at the vicinity of the node. Thus, the analysis 
provides the forces at nodes where displacement is prescribed, and 
displacements where forces are prescribed. If {F„ Uy, is 
prescribed at a surface node for example, one gets {u„ Fy, F^ at 
that node. Even though the solution is unique in principle, errors 
are introduced if the discretization is coarse. With finer 
discretization, the solution converges to the correct one. This 
convergence test is often employed to gage the accuracy of the 
solution. 

In our problem with two cells, we prescribe zero displacement 
boundary conditions at the bottom surface and at the four vertical 
sides of the body (Fig. 2). Thus all the nodes on the bottom and the 
vertical sides are fixed. For simplicity of illustration, consider that 
there are a few nodes on the top free surface outside the cell 
boundary, and a few nodes within (Fig. 2). Our objective is to 
calculate the traction on these nodes. We can experimentally 
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measure displacements {u„ Uj, at all the nodes on the surface. 
They are generated by cell forces, although we do not know the 
precise locations of these forces. We also know that the surface 
outside the cells has no traction, and that each cell or cell cluster 
produces a traction field that is self-equilibrated, i.e., the sum of 
forces applied by the cell or the cell cluster on the substrate is zero. 
Cell traction can be evaluated by prescribing either of the two 
boundary conditions: 

i) Whole-field displacement boundary conditions 

(BCl): [u„ Uy, a J are prescribed at all the nodes on the 
surface, and their traction is analyzed by FEM. 

ii) Mixed boundary condition (BC2): zero traction is 
prescribed at all the nodes outside the cells {F„ =Fy,=F^ = Q), 
and displacements {u„ Uy, are prescribed for nodes within 
the cell boundaries (Fig. 2). 

Remarks. (1) The mix(;(l boundary scheme appli(;s exact 
boundary condition (zero force) at nodes outside the c:ells. Hence 
none of the displacements (m„ Uy, need to be prescribed at these 
nodes. Thus, it is not necessary to measure the displacements of 
the beads outside the cells. Due to the exact boundary conditions 
outside the cells, the traction solution is expected to be more 
accurate. However, errors will be introduced if the cell boundary is 
incorrecdy defined and there are nodes that fall outside the cell 
boundary where cells apply traction. In cases where the cell 
boundaries cannot be identified due to imaging conditions (Fig. 3), 
displacements should be prescribed for regions nearby the cells. 

(2) Displacement and Poisson's ratio: It is shown in the 
supplementary material (Supplementary materials text S3, Fig. 
Sib and c), that if the Poisson's ratio of the gel approaches 0.5, 
then the in-plane displacements, {Uj, u^, on the surface of the gel 
are independent of the out-of-plane component of traction (F^. 
That is, {u„ are determined by {F„ Fy) on the surface. Similarly, 
is determined by F^ on the surface only. Thus, in order to 
evaluate the in-plane traction only, one needs to measure and 
prescribe in-plane displacements only at the surface nodes, and 
prescribe arbitrary boundary condition in z direction (i.e F^ = 0 or 

= 0) at all surface nodes, when Poisson's ratio is close to 0.5. We 
experimentally measured the Poisson's ratio of our gel as 
0.47±0.02 (Fig. S3b, n = 5). In order to estimate the in-plane 
traction only, we have prescribed F^-. = 0 for all nodes within the 
cells in the rest of the paper. This results in an error of less than 
2% in the calculation of in- plane forces F^ and Fy (Supplementary 
materials text S3 and Fig. S3b). If F^ is desired, one needs to 
measure and prescribe {u„ Uy, at the surface nodes. Also, if 
Poisson's ratio is much less than 0.5 (e.g., 0.35), {u„ Uy, must be 
prescribed at the nodes within the cells even when only in-plane 
traction is desired. 

Validation of uniqueness of solution in finite-element 
models 

In this section, we demonstrate computationally that the 
traction solution from finite element simulation is unique as long 
as the fuU 3D boundary conditions are prescribed. We define two 
circular boundaries representing two cells with half-ceU distance 

apart on a soft gel surface. The diameter of each boundary is 
chosen as 20 |J,m, close to real cell size. A three-dimensional finite- 
element (FEM) block model is generated (ANSYS 12.0 Work- 
bench Package) to represent the PA gel substrate [79-98]. The gel 
is presumed linear elastic, isotropic, and homogeneous in their 
mechanical properties for a wide range of deformations [78,99]. 
The Elastic modulus, E, of the gel is IKPa (our experimental value 
is 1.05±0.17 kPa, measured by AFM indentation (n=15; Fig. 



S3a), [99-101]]). The model height is 70 |J,m, same as tiie 
thickness of PA gel used in experiments. We first apply an in-plane 
force field (Fig. 4a) within each boundary, and compute the 
corresponding displacement field, u^, Uy, u^ (Fig. 4b). Second, we 
use the computed u^, Uy and u^ within the cell boundary on the 
surface (Fig. 4c), and zero-traction conditions outside the 
boundaries to calculate the traction within the cells (Fig. 5d). A 
comparison between the prescribed and the calculated forces from 
the two steps shows close quantitative agreement (within 1%) 
(Fig. 4e-f). Note that individual cells or cell clusters generate self- 
equilibrated traction on the substrate. Hence, we use a measure of 
accuracy of the traction solution by defining the error ratio, 

E = ^(-LF^if + Q:Fyif/^J{W.i\f + iWyAf (2) 

where F„ and F^i, are the nodal force components within the 
individual cells, and i = 1 , n, the number of nodes within the cell or 
cell cluster boundary. For exact solution, e = 0. 

Demonstration of 2-cell experiments using mixed- 
boundary condition method 

In this section, we dem()nstrat(; the applicability of the method 
by evaluating the traction induced by two neighboring cells. Here, 
two monkey kidney fibroblasts were plated on PA gel (1 kPa) with 
Poisson's ratio of 0.47 (Fig. 5a). Two different regions (two sets of 

and S„) were selected to prescribe the displacement boundary 
conditions: (1) displacement field underneath the two cells were 
prescribed in the model (the white parts in Fig. 5b), whereas the 
tiaction-free condition was applied outside the cells (the black part 
in Fig. 5b); (2) the displacement field within a region enclosing 
both cells was prescribed (the white part in Fig. 5c), whereas the 
traction-free condition was applied outside this region (the black 
part in Fig. 5c). The out-of-plane force, Fz, was prescribed as zero 
within the cellular regions in (1) and (2). The traction fields were 
calculated for both cases (Fig. 5d,e,g,h), and compared (Fig. 5f and 

i). The RMS =\ r '^'^^^^'y of node-by-node 

V n 

traction difference inside 2-cell region (superscripts indicate 

regions 1 and 2) was 21.7 Pa, which shows close match with only 

5.1% of maximum traction inside the cells (426. 8 Pa). 

Demonstration of whole-field displacement boundary 
conditions method and comparison 

In this section, we compare our mixed-boundary condition 
method with traditional whole-field displacement boundary 
condition method, which requires iterative calculation and has 
been successfully used by Fredberg, et al [49,102]. Briefly, the 
iteration calculation proceeded as follows: (a) we assigned the 
complete 2D DICM (digital image correlation method) displace- 
ment data (ux,Uy) for all nodes of the top surface of the gel (both 
intracellular and extracellular regions; Fig 6a-b). We prescribe 
F^ = 0 within the cluster for both the mixed boundary condition 
and iterative methods, (b) The traction field was solved using 
FEM. Then all the forces in the extracellular region were replaced 
by Fx = Fy = F^ = 0 to satisfy the traction-free condition, while the 
forces in the intracellular region were retained intact, (c) The new 
traction field was used to generate a new displacement field using 
FEM. Thus a new displacement field was computed within the 
intracellular region, (d) The computed intracellular displacement 
field was replaced with the DICM displacement field (u, and Uy), 
while the computed extracellular u^, Uy, and u^ from previous step 
were retained intact, (e) The steps (b), (c), (d) were repeated until 
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Figure 2. Two cells applying contractile forces on 2D elastic substrate. In our finite element scheme, all nodal displacement underneath the 
cells on the top surface of the gel are measured, while for the nodes outside the cells all tractions are assigned zero and thus their displacements are 
not necessary to measure. All nodal displacements at the bottom and side walls of the gel are assigned zero (not shown in the figure). These 
combinations of data inputs constitute the set oi Mixed Boundary Conditions in our FEM simulation. The computed parameters of the model are nodal 
traction underneath the cells as well as displacements of the extracellular nodes (traction-free nodes on the surface). 
doi:10.1371/journal.pcbi.1003631.g002 




Figure 3. The traction for the nodes far from the cells are zero, however errors will be introduced if the cell boundary is poorly 
defined and there are nodes that fall outside the presumed cell boundary where cells may apply traction. In cases where the cell 
boundaries cannot be identified due to imaging conditions, displacements should be prescribed for regions nearby the cells. 
doi:1 0.1 371/journal.pcbi.1 003631 .gOOB 
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Figure 4. Validation of the accuracy and uniqueness of the finite element solution to extract 3D traction force fields, (a) A 

computational model with two regions representing 2 separated cells, each 20 \im in diameter and separated by half-cell distance 10 \im, was 
established. A self-equilibrated force field was applied within each region. The magnitude and directions of forces were indicated by arrows, (b) The 
resultant full displacement field was obtained by ANSYS. (c) The displacement fields underneath each cell were chosen and assigned to the same 
model. The boundary conditions of nodes outside the regions were set traction-free, (d) A new force field was obtained using the above mixed- 
boundary condition. The magnitude and directions of nodal forces were shown by arrows, (e-f) The node-by-node difference between initially 
applied forces and retrieved forces (in x and z direction, respectively) are shown. The difference is <10^^ nN (within 1%). 
doi:1 0.1 371/journal.pcbi.1 003631 .g004 



the solution converged, i.e., the difference between the root mean 
square (RMS) of surface nodal forces in two consecutive cycles 
became less than 5% (Fig. 6c-e). 

Our computational results showed that the solutions from 
mixed-boundary and iterative methods converge (Fig. 6c-e). We 
found, the difference between the root mean square (RMS) value 
of traction from the two methods was 1.6x10"' kPa (Fig. 6f), less 
than 3.8% of the maximum computed cell traction. The difference 
between the RMS of the nodal forces was 0.2 nN, which is 0.25% 
of the maximum nodal force at cell cluster - substrate interface 
(Fig. 6g). The distribution of traction |t| and forces |F| at nodes 
(Fig. 6h-6j) shows good agreement between the two methods. We 
used E (Eqn 2) as a measure of accuracy of the traction solution. 

Mesh size effect 

In FEM, convergence test is required to determine the optimal 
mesh size needed to obtain the accurate solution. Three mesh 
sizes, 3.23 |J.m, 4.84 |j,m, and 6.45 |J,m were tested, as shown in 
Fig. 7a-c, and used to calculate the traction field of the same cell 
cluster by mixed-boundary condition method. The distribution of 
nodal traction and forces showed minor difference between the 



three mesh sizes (Fig. 7a-c and 7e). The values of s were 4.74%, 
6.69%, and 6.12% for mesh size of 3.23 nm, 4.84 |xm, and 
6.45 |j.m respectively (Fig. 7d). Therefore, in the following 
computations, mesh size of 4.84 |j.m was used for analysis. The 
upper limit of mesh size is dependent on the specific cell size and 
the gradient of the traction field produced by the cell. A starting 
point on mesh size can be <20% of cells size. 

Traction calculation for multiple cell clusters 

A key attribute of the present method is the computation of 
traction fields generated by multiple cell clusters interacting with 
each other. Each cluster may consist of multiple cells, and the 
cluster size might be similar to or larger than the thickness of the 
soft substrate. Hence the effect of the glass-gel interface needs to be 
considered, and the gel may not be treated as half space. In the 
following, we study several cell clusters (Figs. 8-10) and outline the 
main biological findings. The mixed-boundary condition method 
was used to compute the traction fields. 

Cells in clusters exert cell-cell forces. Fig. 8a shows a 
cluster of colon cancer cells (HCT-8) on 2 kPa substrate. These 
cancer cells are epithelial in nature, and have low metastatic 
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Figure 5. Verification of the uniqueness of solution of the traction field computed from the experimental displacement field, (a) 

Phase-contrast pictures of 2 spatially isolated MKF cells on 1 kPa PA gels, cultured after 1 day. Scale bar: 15 |im. (b) The displacement fields 
underneath each cell were chosen for computation, (c) A larger area enclosing both the cells and neighboring area was chosen where displacement 
field was prescribed. In both cases, the nodes outside the selected regions were set traction-free, (d-e) The traction field computed by above 2 cases 
were visualized and compared by 2D contour plots (d-e) and 3D bar representation (g-h). Also, the node-by-node difference of traction fields 
computed using 2 selected schemes was illustrated by both 2D contour plot (f) and 3D bar representation (i). Dashed lines in orange outline the cells 
boundaries. 

doi:1 0.1 371 /journal.pcbi.1 003631 .g005 



potential. Figs. 8b and c sliow the traction generated by the cluster 
and the corresponding nodal forces on the substrate. Here the grid 
size in the analysis is about 5 |J,m. W e note that the cells within the 
cluster do not form individual force dipoles. The entire cluster 
behaves as a single contractile unit, and generates forces along the 
periphery. The concave peripheries generate larger inward forces. 
These forces are balanced by far away opposing forces. Thus the 
cells behave as generators and transmitters of forces from one edge 
to the other of the cluster. The mechanics of this force 
transmission can be visualized from a free body diagram 
(Fig. 8d). Here the cluster exerts contractile forces on the substrate 
through the acUiesion sites of the outer cells. The cells inside the 
cluster contribute and transmit the force possibly through cadherin 
junctions and force-bearing cytoskeleton. Thus the cell-substrate 
traction is transduced to cell-cell contractile forces. As if, the 
peripheral cells puU the interior cells outward. This is consistent 
with the observation that advancing edge of a cell cluster puUs the 
interior cells [49,102,103,104,105,106,107]. Previously we report- 
ed that HCT-8 clusters are sensitive to their mechanical 
microenvironment and display a metastasis-like phenotype 



transition (MLP) when cultured on appropriately soft substrates 
[7,8,28,29,30,31]. This MLP transition always initiates from the 
periphery of clusters. It remains to be seen whether the difference 
in forcing on the peripheral cells compared to those in the interior 
contributes to the MLP transition. 

Cell clusters may generate traction interior to the 
periphery. Fig. 9a shows two clusters of pre-MLP HCT-8 cells 
close to each other, cultured on 2 kPa substrate for 5 days. The 
traction and the force maps (Fig. 9b, c) show that there are two 
regions well within the boundaries where the traction is high, 
unlike the cluster of Fig. 8 where the traction was mostiy 
peripheral. Here the interior forces are balanced locally, i.e., these 
forces form local dipoles leaving the rest of the cluster nearly 
traction free. This may explain the spherical morphology of the 
clusters, which minimizes their surface areas. The cells of the 
clusters might be under compression due to line tension of the 
peripheral cells. 

Traction domains of cell clusters are dynamic. The cell 

clusters of Fig. 9a merge on the 6th day of culture, as shown in 
Fig. 9d. Soon after merging, the traction and force pattern changes 
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Figure 6. Comparison of mixed-boundary condition metKiod and full-field displacement boundary condition method, (a) Phase 
contrast picture of a single cell cluster to be studied. Scale bar: 50 ^im. (b) The displacement field generated by cell cluster on the top surface of 
substrate, (c-e) The traction field calculated by mixed-boundary method, and full-field displacement boundary method (with iterative calculation 1 
time and 2 times, respectively), were shown respectively. The difference of RMS of the traction between mixed-boundary method and full-field 
displacement boundary method with 1 time iteration was 1.6x10^' kPa, less than 3.8% of the maximum computed traction at cell cluster and 
substrate interface. The difference of RMS of their nodal force was 0.2 nN, which was 0.25% of the maximum nodal force at cell cluster and substrate 
interface. Dashed lines in orange outline the cells boundaries, (f-g) Histograms of nodal traction and force obtained by the two methods 
demonstrated good agreement between each other, (h) Sum of net forces and absolute forces calculated by the above three conditions. The force 
equilibrium was best satisfied in mixed boundary condition method, which is 6.69% of total force, (i) Sum of surface nodal force distribution 
calculated by above three conditions. The RMS results of nodal force calculated by mixed BC method and 1-time iteration method agreed within 
4.96%. (j) Sum of surface nodal traction distribution calculated by the above three conditions. The RMS results of nodal force calculated by mixed BC 
method and 1-time iteration method agreed within 9.27%. 
doi:10.1371/journal.pcbi.1003631.g006 
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Figure 7. (a-c) The convergence test was performed to determine the maximum fine mesh size needed to obtain the accurate 
solution. The mesh sets with different Ax and Ay (Ax = Ay = 3.23 |im, 4.84 \im, and 6.45 \im, respectively) were tested respectively. The traction 
distribution map and traction magnitude histograms from three mesh-size displayed uniform feature patterns. Dashed lines in orange outline the 
cells boundaries, (d) All three cases showed sum ratio of net forces within 7%, satisfying the force equilibrium requirement, (e) The root mean square 
(RMS) difference of traction between 3.23 and 4.84 ^im meshes was about 64.06 Pa (1.28% of maximum computed traction), and the difference 
between 4.84 and 6.45 ^im mesh sizes was about 192.7 Pa (3.86% of maximum traction). The comparison indicates that when mesh size is reduced to 
4.84 |j,m or below, the traction output starts to show minimum variation. 
doi:1 0.1 371/journai.pcbi.1 003631 .g007 



dramatically (Fig. 9e-f). The net force increases by about 20 folds, 
although the direction of the net force does not change. Many 
more cells in both the clusters now participate in generating 
traction. The new traction regions are also mostly interior to the 
periphery, and the merged cluster takes a smooth elliptical shape. 
The merger between two cell clusters mimics that between two 
droplets, as they both tend to minimize the surface energy. It is 
known that cells may interact with each other through substrate 
strain fields [4,11,108,109,110]. In case of the two neighboring 
clusters, the displacement fields were localized well within the 
clusters. It is thus unlikely that their merger was induced by strain 
fields. 

Evidence of cell-cell compression. It is generally under- 
stood that cells generate contractile forces produced by actomyosin 
machinery [11,76,111,112]. However, in a 2D cluster, cells may 
be subjected to compression as shown in Fig. 10. Here, monkey 
kidney fibroblasts (MKF) form several large and small clusters on 
1 kPa substrate. Each cluster is sufficiently far away from the 



others so that there is no mechanical coupling between them. The 
displacement field between them is negligible (Fig. 10b). The 
traction within each cluster is shown in Fig. 10c. Unlike the 
previous two examples, here many more cells in the clusters 
participate in traction generation. Fig. lOd shows the nodal forces 
of the larger cluster. Here, several regions generate dipole type 
forces within the cluster. However, there are interior boundaries 
where opposing forces appear on the substrate, i.e., cells "push" 
against each other. This can be explained by the schematic of 
Fig. lOe where neighboring cells have adhesion sites with the 
substrate. Due to the low stiffness of the substrate, the cells have 
less likelihood of spreading or wetting the substrate, though they 
may adhere to the substrate due to the fibronectin functionaliza- 
tion. Now, if growth occurs in any of the cells next to a neighbor, it 
would push out the neighbor generating an outward force on the 
substrate. This results in a compressive stress between the cells. 
There are three regions of such cell-cell compression in the cluster 
of Fig. lOd (enclosed by dashed lines). 
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Figure 8. Traction and force maps of single Kiuman colon cancer cell (HCT-8) cluster. The cluster behaves as a single contractile unit, (a) - 
(c): Phase-contrast image, traction and nodal force map of a well-spread pre-IVILP HCT-8 cancer cell cluster. The cells were cultured on 2 kPa hydrogel 
substrates. The distance between the nodes is about 5 |im. Scale bar: 60 |im. Colors of contour represent the magnitude of traction stress. Vectors 
indicate the direction of traction force at each node and arrow lengths represent the magnitude of node force. Dashed lines in orange outline the 
cluster boundary, (d) A free body diagram visualizes the mechanics of this long-distance force transmission. The cell cluster exerted contractile force 
on the substrate through the adhesion sites of the outer cells. The inner cells transmitted the force possibly through cell-cell junctions and cortical 
actin. 

doi:10.1371/journal.pcbi.1003631.g008 



Discussion 

The majority of fundamental physiological processes in tissue 
development, health, and disease are coordinated by the collective 
activities of multiple cells [60,62,76,102], rather than single 
cells[10,103]. To understand how mechanical traction applied 
by neighboring cell cluster groups could specify or mediate the 
tissue functionalities [7,8,11,75,104,105,106], robust cellular 
traction evaluation method is indispensable. In the present study, 
we developed a finite element element-based traction force 
microscopy (TFM) to accurately compute and visualize the 
traction maps resulting from multiple cell clusters. The uniqueness, 
convergence, and correctness of traction solutions are substanti- 
ated. We showed that as the gel Poisson's ratio >0.4, the in-plane 



traction can be obtained with minimal error from the in-plane 
displacement field alone. For Poisson's ratio <0.4, both in and out of 
plane traction depend on both in and out of plane displacement 
boundary conditions, and it is essential to measure these displacements 
to compute any of the traction components. The method presented is 
appficable to substrates widi any value of the Poisson's ratio. It 
calculates the fiill 3D traction field given the 3D displacement 
boundary condition within cells or cell clusters. Moreover, unlike the 
classical TFM metiiods that are based on Boussinesq solutions 
[39,40,48,49], the FEM takes into account the effect of substrate 
thickness and nearby environment. It is now known that cells can 
sense the substrate depth within the cellular length scales by showing 
distinct morphological variation on the gel substrate with same Elastic 
modulus but with varying thickness [22, 107]. 
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Figure 9. Traction maps of two neighboring human colon cancer cell (HCT-8) clusters. Their interior traction domains are dynamic, (a) - (c): 
Phase-contrast image, traction and nodal force maps of two independent cancer cell clusters cultured on 2 kPa flexible hydrogel. Cells were on 
culture day 5. Each cluster generated high traction well within the periphery, leaving the periphery almost traction-free, (d) - (f): Phase-contrast 
image, traction stress and nodal force maps of the merged pre-IVILP HCT-8 cancer cell cluster after 24 hours (6* culture day). Following merging, 
many more cells in both the clusters participated in generating traction, and the net force increased by about 20 folds, although the direction of the 
net force did not change. Scale bar: 40 |im. Dashed lines in orange outline the cluster boundaries. 
doi:1 0.1 371 /journal.pcbi.1 003631 .g009 



We applied the method to compute the traction generated by 
multiple cell clusters. Some of the clusters were more than 100 |j.m 
in size consisting of many cells, while others were in close 
proximity to each other. The computational scheme presented 
here is ideal for studying such clusters, since the domain of traction 
field is much larger than the thickness of the gel, and one needs to 
account for the finite thickness of the substrate. A few interesting 
biological insights emerge from these analyses. First, the cluster 
may behave as a single contractile unit where the peripheral cells 
serve as anchorage sites. Force is transmitted between distant 
peripheries by the cells inside the cluster. Thus the cells are 
subjected to tensile intercellular forces, as if the peripheral cells are 
pulling the interior cells outward. It needs to be seen whether there 
are specific cells within the cluster that generate the force, or all 
the cells behave as contractile actuators. In any case, the cells 
probably use cell-cell junctions and cytoskeleton to transmit the 
force through the cluster. 

We also found instances where traction is limited to small 
regions well within the clusters. These regions can have locally 
balanced traction (forming dipoles), leaving the rest of the clusters 
nearly traction free and weakly adhered to the substrate. These 
clusters are spherical in morphology, as expected. The traction 
free regions tend to minimize the surface area by being circular, 
just as a free-standing cell cluster takes a spherical shape. It is 
plausible that the cells within the circular clusters are under 
compression due to the surface tension of the peripheral cells. In 
any case, the interior traction maps can be highly dynamic. When 
cell clusters merge, the traction map can change their orientations, 
and the net force can increase by an order of magnitude over short 
times. 

It is known that cells generate contractile forces. Thus, it is 
expected that the cells in a 2D cluster will be under intercellular 



tension. We found evidence to the contrary. If the cells are on soft 
substrates where they do not spread much, but they adhere to the 
substrate, then some of the cells in the cluster may be subjected to 
compression. We found regions within such clusters where the 
neighboring cells apply repulsive forces on the substrate, i.e., the 
cells are pushing against each other while being adhered to the 
substrate. One possible explanation might be that the neighboring 
cells are growing, but their adhesion sites are stationary. 

In conclusion, we developed a robust FEM-based cell traction 
force microscopy technique to estimate the traction forces 
produced by multiple cells and clusters. The utility of the 
technique is exemplified by computing the traction force fields 
generated by multiple monkey kidney fibroblast (MKF) and pre- 
MLP human colon cell (HCT-8) clusters in close proximity. The 
developed technique is user-friendly and computationally inex- 
pensive. Our FEM-based traction force microscopy provides a 
powerful tool to probe multi-cell questions involving assembly/ 
disassembly dynamics of cell ensembles, tissue network formation, 
and wound healing. Future work is needed to determine the 
subcellular processes involved in mechano-sensing and regxilation, 
and their respective timescales. 

Materials and Methods 

PA gel substrate preparation and ECM conjugation 

Polyacrylamide (PA) gel substrates with 1 kPa stiffness used in 
present study were made by mixing 12.83% (v/v) of acrylamide 
(Sigma-Aldrich, Inc.), 1.54% (v/v) of N, N-methylene-bisacryla- 
mide (Sigma-Aldrich, Inc.), 2% (v/v) of 1 |J,m diameter fluorescent 
micro-beads (Invitrogen, Inc.) and 10 mM Hepes (Gibco., Inc.) 
[7,11]. Solution was vortexed thoroughly for 5 min to obtain 
uniform distribution of beads. TEMED and ammonium persulfate 



PLOS Computational Biology | www.ploscompbiol.org 



11 



June 2014 | Volume 10 | Issue 6 | e1003631 



The Study of Forces within Multi-cellular System 




Traction r\ 





Figure 1 0. Evidence of cell-cell compression in monkey kidney fibroblast (MKF) cluster, (a) Several MKF cell clusters on 1 kPa PA gel. (b)- 
(c) The displacement and traction field produced by the clusters on the top surface of the substrate. The traction by the small clusters is negligible 
compared to those generated by the larger ones. Dashed lines in orange outline the cluster boundaries, (d) Nodal forces computed for the largest 
cluster. The finite element grid size is about 5 \irr\. There are regions in the cluster, shown by dashed lines, where repulsive forces appear on the 
substrate, i.e., cells "push" against each other, (e) To explain the cell-cell compression, a free body diagram is shown to reveal the intercellular force 
and cell-substrate traction force of 2 neighboring cells on the substrate. As the substrate is soft, the cells have less likelihood of spreading or wetting 
the substrate, but can adhere to the substrate due to the fibronectin functionalization. As cell proliferation and growth occur within the cluster, the 
cells push against their neighbors, generating an outward force on the substrate. Scale bar: 50 |im. 
doi:10.1371/journal.pcbi.1003631.g010 



(Fisher Scientific, Inc.) were used to initiate PA gel crosslinking. 
Chemical modification of glass slides and preparation of PA gels 
were carried out following the procedures described previously 
[50,80,81,82,83,84,85]. Briefly, a circular glass coverslip (Fisher 
Scientific, Inc.) of 1.2 cm in diameter was placed on an acrylamide 
solution drop on activated coverslip and placed on the bottom of a 
petri dish. Capillarity spreads the drop and fills the space between 
the circular coverslip and the activated coverslip. The gel was 
cured at room temperature and reached to the stabilized thickness 
of 70 |J.m [82,85,86]. The circular glass coverslip was peeled off 
from the gel that remained on the activated cover slip. The 
surfaces of the air dried PA gels were activated by incubating in 
97% hydrazine hydrate (Acros Organics.) for 12 h followed by a 
complete rinsing with DI water and 30 minutes incubation along 
with gentle shaking in 5% acetic acid (Avantor Performance 
Material, Inc.) [7,8,11,13,81]. Solution of human fibronectin 
(25 |a,g/ml, BD Biosciences) was prepared by dissolving in 
phosphate buffer saline (PBS) and the carbohydrate groups of 
fibronectin were oxidized by sodium periodate (Sigma-Aldrich, 
Inc.). To minimize the displacement noise and rigid body motion 
during imaging, the glass slides was firmly adhered to the bottom 
of 30 mm petri dish using adhesive glue (Henker Consumer 
Adhesive, Inc.). Full experiment procedures and sample charac- 
terization are provided in Supporting Materials Text S4-S9 and 
Figures S1-S5. 



Supporting Information 

Figure SI Gonfocal microscopy images of monkey 
kidney fibroblasts (MKF) cells on gel substrate with 
immunofluorescent stained F-actin cytoskeleton (green) 
and focal adhesion protein, Vinculin (red). The x-y plane 
shows the horizontal view of spread MKF cells. The z-y and z-x 
cross-sectional planes, A and B, show the vertical structures of 
spread MKF. They display that the height-to-length ratio of 
spread MKF cells is in the range of 1 /40~ 1 /50. Vinculin staining 
(red) indicates the basal surface of MKF cells. This low height-to- 
length ratio implies that the cells exert their traction forces mosdy 
along the x-y plane through their contractile filaments. The 
cartoon of spread cells on top right of (a) shows that the angle (p 
between contractile cytoskeleton (green) and substrate is within the 
range of 1~10". (b) To estimate the error due to out-of-plane 
forces on the evaluation of in-plane traction, a general 3D force- 
displacement model for the cell is developed. In the model, the cell 
applies both in-plane and out-of-plane forces on the substrate, 
and P, with corresponding deformation u and w. (c) The error 
index plotted for all three cases, P — 0, w = 0, and general loading 
P = kQ^v.s Poisson ratio v. For P—0, there is no error in planar 
force calculation for all position x and displacement u. For other 
cases, however, there are deviations due to the presence of out-of- 
plane force, P, at different boundary conditions. It is evident from 
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the plot that for small values of Poisson's ratio, the z-component of 
deformation w will influence the in-plane force (),and thus create 
varying results depending on loading modes and the value of 
Poisson's ratio. Therefore, excluding out-of-plane deformation w 
will introduce error in calculating the in-plane force, Q. However, 
as Poisson's ratio approaches Vi, most of the discrepancies in 
planar force calculations becomes negligible, and all set of curves 
converge to a unified value corresponding to P=0, regardless of 
value and direction of the out-of-plane force P and deformation w. 
(TIF) 

Figure S2 A representative elastic body subject to the 
most general form of mixed boutidary condition. 

Displacement field and traction field are given on separate 
surfaces and S„ respectively. The body has total surface 
S = S„+S„ and total volume V. The general state of stress tensor ay 
and rc-spective displacement vector Ui are shown at an arbitrary 
point P within the body. Cauchy traction vector f,- applies on an 
arbitrary, infinitesimal surface denoted by the unit normal vector 



Figure S3 Measurement of PA gels' Young's modulus 
and Poisson's ratio, (a) The PA gel stiffness was measured by 
AFM as 1.05 ±0.1 7 kPa (n = 15), and fitted by Hertz's indentation 
theory, (b) Uni-axial tension experiments were carried out to 
stretch PA gel samples with dimension 2.2 cmx5.0 cmx4.0 mm 
under aqueous condition. The lateral and axial strains were 
recorded progressively and fitted into a linear plot to obtain the 
Poisson's ratio. The Poisson's ratio was determined as 0.47 ±0.02 
(n = 5) and appeared to be independent of gel bulk stiffness. Two 
representative examples are shown. 
(TIF) 

Figure S4 Contour plots show the displacement field 
produced by the MKF cell obtained by a commercially 
available DIC software VIC-2D (a) and by the open 
source MATLAB DIC program (b), respectively, (c) The 

node-by-node displacement difference plot shows that the two 
DICM methods give quantitatively similar displacement data. 
(TIF) 

Figure S5 (a) A Tungsten probe with known stiffiiess of 
10.74 nN/|tm (calibrated with weight) was vertically 
held by a high-resolution x-y-z piezo-stage to apply 
horizontal force on the flexible hydrogel surface, (b) The 

deflections of probe tip with respect to reference base, as well as 
the resultant displacement fields of beads on gel's top surface, were 
recorded. The displacement fields were assigned to FEM model to 
compute the resulting force. The double-headed arrows indicated 
the gap between micro-needle and reference base. Multiplying this 
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